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Abstract 


A computational plasma aerodynamics model is 
developed to study the performance of a laser propelled 
Lightcraft. The computational methodology is based on 
a time-accurate, three-dimensional, finite-difference, 
chemically reacting, unstructured grid, pressure-based 
formulation. The underlying physics are added and 
tested systematically using a building-block approach. 
The physics modeled include non-equilibrium 
thermodynamics, non-equilibrium air-plasma finite-rate 
kinetics, specular ray tracing, laser beam energy 
absorption and retraction by plasma, non-equilibrium 
plasma radiation, and plasma resonance. A series ol 
transient computations are performed at several laser 
pulse energy levels and the simulated physics are 
discussed and compared with those of tests and 
literatures. The predicted coupling coefficients for the 
Lightcraft compared reasonably well with those ot tests 
conducted on a pendulum apparatus. 

Nomenclature 

Aj reaction rate constant of the j-th reaction 

Afj.Am matrix coefficients of transport equations 
c speed of light in a vacuum. 
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vibrational energy 
equilibrium vibrational energy 

flux vector 

Gibbs’ free energy of species i 
total enthalpy 

intensity or number of quantum levels 
emission coefficient 

transfer constant between quantum levels 

Boltzmann’s constant 

Coulomb constant 

thermal conductivity 

molecular weight 

electron mass 

total quantum level population 
quantum level 
refraction index 
unit normal vector 
electron number density 
Prantdl number 
gas static pressure 
partial electron temperature 
collision energy transfer 
net radiative heat flux 
electron charge 
gas constant 
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r 

position vector 

S 

source terms 

T 

gas temperature, K 

T e 

electron temperature, K 

T ( 

translational temperature, K 

T v 

vibrational temperature, K 

t 

time, s 

lit 

velocity components 

V 

velocity magnitude 

Xt 

Cartesian coordinates 

Y 

i 

chemical species mass fraction 

Peek 

Symbols 

eo 

permittivity 

0 

dependent variable 

r 

control surface 

Kl 

laser absorption coefficient. 

Xe 

electron species thermal conductivity 

Xi 

laser wavelength 

P 

viscosity 

Vi/ 

stoichiometric coefficient 

e 

refractive angle 

p 

density 

r it 

vibrational-translational relaxation time scale 

TV 

shear stress tensor 

Q 

control volume domain 


chemical reaction source term 

top 

plasma frequency 


Introduction 

Beamed energy propulsion was tirst promoted by 
Kantrowitz 1 . Since then, a propulsion system supported 
by a laser-sustained plasma has been the subject of many 
studies.' 6 The main advantage gained by laser 
propulsion is the low-weight system derived from 
decoupling the energy source from the vehicle, and high 
specific impulse resulting in low tuel consumption. The 
first ground and flight tests of a ground-based laser 
propelled vehicle were reported in 1998, 7 in which active 
tracking and beam control were demonstrated to 122 m 
on a horizontal wire and spin-stabilized tree flights in the 
laboratory were accomplished to altitude of 4 m. Six 
months later, the spin-stabilized vertical free flights 
outdoors reached 30 m. These results were 
accomplished with axisymmetric vehicles which are of a 
special design in which a nosecone shaped forebody, an 
annular shroud, and a parabolic afterbody are the only 
major components. A vehicle of such design is hereby 
named as the “Laser Lightcraft’’ in this study. The layout 
of a computational grid for a Laser Lightcraft is shown in 
Fig. 1. The parabolic afterbody surface serves both as an 
aerospike nozzle and also the main receiving optic, 
whereas the annular shroud surrounds the ring focus of 
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the parabolic optic. The Laser Lightcraft uses a Rocket- 
Based Combined Cycle (RBCC) to operate in both the 
rocket and air-breathing modes. When operating in an 
air-breathing mode, the specific impulse is infinity since 
no fuel is consumed. These successful tests 
demonstrated the concept and the potential feasibility of 
launching small payload with a Laser Lightcraft. 

When the focused beam energy strikes the shroud of a 
Laser Lightcraft, free electrons form to invoke Inverse 
Bremsstrahlung (IB) and the optical breakdown ensues. 
In order to achieve maximum intensity and to avoid 
plasma resonance, repetitive short pulses are used. As 
such, repetitive Pulsed Detonation Waves (PDW) are 
generated to propel the vehicle, hence the Laser 
Lightcraft engine is also known as a PDW engine. It can 
be seen that the underlying physics involved such as IB 
and PDW are so complicated that a simple system model 
is not capable of describing the phenomena, let alone 
predicts the thrust performance and the associated flow 
and thermal environments necessary for designing for the 
structure integrity. A detailed computational plasma 
aerodynamics model is therefore needed for the design 
and scaling of future Laser Lightcrafts. 

From May to July, 1999, performance data in the form 
of coupling coefficients were taken in a series of indoor 
tests for the Model #200 series vehicles. These test were 
performed under a cooperative effort between Rensselaer 
Polytechnic Institute, the Air Force Research 
Laboratory's Propulsion Directorate and NASA’s 
Marshall Space Flight Center (MSFC), to enhance Laser 
Lightcraft propulsion research. These data are the first of 
a kind benchmark quality data for computational model 
development. Armed with these data tor comparison 
purposes, a computational plasma aerodynamics 
methodology is developed as a first attempt to predicting 
the thrust performance of a PDW engine propelled Laser 
Lightcraft. 

rnmpntntional Plasma Aerod ynamics Development 

Several attempts have been made to analytically or 
computationally model the physics inside continuous 
wave (CW) devices. For example, an early effort by 
Raizer 9 assumed constant-pressure one-dimensional (ID) 
flow' in air with heat addition by a laser. Thermal 
conduction was considered dominant and the radiation 
loss from the plasma was ignored. Kemp and Root 
later extended this ID analysis to hydrogen and to 
include thermal radiation. Molvik, Choi, and Merkle 
extended this problem to a two-dimensional (-D) 
structured-grid hydrogen flow and implemented real ray 
tracing, but ignored thermal radiation. Jeng and Keefer 
did similar analysis and added thermal radiation. 
Myrabo, Raizer, and Surzhikov 6 later simulated the 
radiative-gasdynamics processes in optical discharges 
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maintained in a subsonic De Laval nozzle. In the work 
described in this paper, a multidimensional unstructured* 
grid computational plasma aerodynamics methodology is 
developed for both CW and repetitively pulsed wave 
devices, with emphases on nonequilibrium effects and 
laser-induced physics such as laser-plasma interactions. 
Realistic laser absorption coefficients are used. The 
Governing equations will be described first, with the 
auxiliary equations that compute the nonequilibrium and 
laser induced physics to follow. 


-^-j Q p0dQ+^ F ndV=j Q S 0 dQ 
at 

where Q is the domain of interest and V denotes the 

surrounding surfaces; n is a unit normal vector of T in 

the outward direction. The flux function ^ contains the 
inviscid and the viscous flux vectors. 


ThP Governing Equations 

To properly describe the plasma aerodynamics 
involved in laser propelled propulsion, the time-varying 
transport equations of continuity, species continuity, 
momentum, global energy (total enthalpy) and electron 
energy are formulated and written in a Cartesian tensor 
form: 


Let’s consider a control volume interface, e, between 
control volumes E and P, with^ normal vector n. And, 

0E -$P ~ # fe ~ 'p), where is interpolated from 

the neighbor cells E and P. For the face e between 
control volumes P and E , the diffusive flux can be 
approximated as: 
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where the shear stress n y can be expressed as: 




f d a, d it 2d u 
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dx, dx. 3 d y 


The beam energy absorption, beam energy transfer, and 
plasma radiation are invoked through the source terms of 
the global energy and the electron energy equations, 
whereas the plasma species are generated through the 
source terms of the species continuity equations. 


{V<p»n) e 


' 4>p +V0 t . 

in ?E ?p . \ 

r p 

r E- r P\ 


The finite volume formulation of flux integral can be 
evaluated by the summation of the flux vectors over each 
face. 


|> F ndY = X C; AT j 

Jr j = tin 

where k(i) is a list of faces of cell i. F, y represents 
convection and diffusion fluxes through the interface 
between cell i and j, and is the cell-face area. 

The convective flux is evaluated through the upwind- 
cell quantity by a linear reconstruction procedure to 
achieve second order accuracy: 

<P e = (p u + • ( r e “ r u ) where the subscript u 

represents the upwind cell and y/^ is a limiter used to 
ensure that the reconstruction does not introduce local 
extrema. The limiter proposed by Barth is used here. 
Defining <Wv = wtxWu j > , ^min = ftun{0 u ^ j ) (and 

assuming $ e is computed with y/ e - 1) the scalar 
associated with the gradient at cell u due to edge ^ is: 


Numerical Scheme 

The cell-centered scheme is employed to allow the 
volume surfaces to be represented by the grid cell 
surfaces. The transport equations can be written in 
integral form as: 



if o 0 > 0 

if 0 ] - 0 < 0 
if 0-0 = 0 
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Solution Pr ocedures 

A general implicit discretized time-marching scheme 
for the transport equations can be written as below. 


equation, the source term is linearized to result in an 
explicit term and an implicit term. This treatment is 
important for an unconditionally stable solution of the 
equation. That is. 


A / 




i n+i 

4>r 


NB 

= 1 4 , 

m~\ 


(PgpT 

At 


+ 5* 


where NB means the neighbor cells of cell P. The high 
order differencing term and cross diffusion term are 
treated using known quantities and retained in the source 
term and updated explicitly. The second term of the 
right-hand side (RHS) is a perturbation that comes from 
the RHS. A predictor and corrector solution algorithm is 
employed to provide coupling of the governing 
equations. The discretized finite-volume equations form 
a set of linear algebraic equations, which are non- 
symmetric matrix system with arbitrary sparsity pattern^ 
The preconditioned Bi-CGSTAB and GMRES(w) 
matrix solvers are used to efficiently solve the linear 
algebraic equations. 


Auxiliary F qnations 

Thermal N o n-Fqnilihriiim F.nerev Fol ia tio ns 

For high temperature flows, thermal non-equilibrium 
state may be important. In Landau and Teller s 
derivation, 14 a master equation is employed to describe 
the evolution of the population of quantum level N„ This 
master equation is written as: 

‘Jlh- - ,vV K^.N - N I K, mi N, :i - 0.1.2 

Ji ,“o i =o 


4 q h) 

Explicit source term = P tLJ 
P 

Implicit source term = TLT 

The vibrational temperature is used to influence the 
reaction rates of chemical reactions by assuming that the 
rate coefficients for dissociation are functions of the 
geometrical mean temperature between T and T v * 

Mnn-Fqiiiiihrium Ai r Ch emistry 

A general system of chemical reactions can be written 
in terms of its stoichiometric coefficients (v|j and vjj ) 
and the i h chemical species name (Af/) of the f reaction 
as 


i i 


The net rate of change in the molar concentration of 
species / due to reactions j, and the species 

production rate can be written as: 




py , : 


Results from the quantum mechanical solution of the 
harmonic oscillator are used to relate the various 
quantum transition rates to one another, and then the 
master equation may be summed over all quantum states 
to obtain the Landau-Teller equation: 


Dpe v 

Dt 


d 



dx t \ 

k V J 

d - X > > 

+ H 

r LT 


The forward (Arrhenius law) and backward reaction 
rates for each reaction is given by: 


Kf } = Aj'Expi—^;) 


Kb, 


K/, 

Kt 

J 


where Kej is the equilibrium coefficient 


An empirical expression (given below) is used to 
model the Landau-Teller relaxation time scale. The 
Landau-Teller expression for the vibrational relaxation 
time is calculated using the 6-coefficient specified. 

_ exp[(fej Iff' j 
T/ r p{\ - exp(- b b !T)) 

where b,=7.1E-4. b : =0, b.i=1.9LE6. b 4 =0.3333, b 5 =0, 
b 0= 0 are used here. To solve this vibrational energy 


Ke, =(RT) 


Exp j - 


n 

V 

1=1 


(v. 


V v'*i 


RT 


i j 

A point-implicit (operator splitting) method is 
employed to solve the chemistry system. For the 
breakdown of air. Park's multitemperature air chemistry 
is baselined in this study. This mechanism composes of 
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the dissociation, NO exchange, associative ionization 
charge exchange, electron-impact ionization, and 
radiative recombination reactions. With this mechanism, 
electrons are produced first by the associative ionization 
process and more electrons are produced by the electron- 
impact reactions. Since the number of electrons doubles 
in each such event, electron density increases 
exponentially, in the form of an avalanche. This 
mechanism thus provides the initial electron density tor 
ignition and produces the avalanche of electrons 
necessary for the subsequent optical breakdown. It can 
be seen that the IB process is embedded in the 
dissociation, associative ionization, and electron-impact 
ionization reactions. The reaction rates were validated 
with experiments of postshock temperatures ranging from 
20,000 - 60,000 K which are inline with the maximum 
computed breakdown temperatures. 
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the energy absorbed by the plasma reaches 15% of that 
delivered. Spark times computed based on that criterion 
range from 0.4 ps to 1.2 ns. As expected, lower laser 
pulse energies require longer spark times for ignition. It 
is also reasonable to assume that the plasma initiates in a 
(spark) region enclosed by the two outer laser rays and an 
arc with an origin at the focal point (in a 2D sense), an 
the (spark) radius of that region is determined empirically 
as 1 mm. In actuality, the spark region is a torus 
venerated by a pie-shaped cross-section rotated about the 
vehicle axis. Finally, the spark power is also empirically 
determined as 15% of that of the incident laser - its 
theoretical maximum. 

Laser Radiation , 1 , 

Geometric optics is used to simulate the local 

intensity of the laser beam, which is split into a number 
of individual rays. In the presence of absorption, the local 
intensity of each ray follows the Beer’s law: 


j dh 

ds, 


~K>1, 


The electron is the only plasma species that a ^orfc’ s 
the laser energy. The absorption coefficient of the CO; 
laser radiation, corrected for stimulated emission m the 
single ionization range, is approximated by the formula: 


K ico , 


5.72p/ln[27<T/10 4 ) 4/3 "' l ' 3 i 


(7710 4 ) 7 ' 2 


plasma Initiation 

Understanding the mechanisms responsible for plasma 
initiation (ignition) has been the subject of ongoing 
research. A spark gap, an extremely intense pulse of 
laser energy striking the focal point, seeded molecules or 
particles, or a retractable tungsten target placed at the 
focal point 16 have been used as sources of tree electrons 
for initiation. 2 It has also been shown that plasma can be 
ignited quite easily ott metal surfaces. Likewise in 
numerical modeling, initial free electrons are required tor 
plasma ignition. This means determining an initial 
electron density or a threshold breakdown intensity. For 
example, Ergun 16 described an initial electron density 
that is necessary for modeling the ignition of a hydrogen 
plasma initiated with a retractable tungsten wire. In this 
study a threshold breakdown intensity approach is 
devised in anticipation that a different threshold 
breakdown intensity may be encountered tor ditfenng 
surface material, target approach, and laser povver. The 
advantage of this approach is that the threshold 
breakdown intensity can either be calibrated or measure . 

For example, the effect of seeds on laser breakdown 
intensity has been reported. 2 In this study, this approach 
is accomplished by performing parametric studies on 

plasma initiation with spark time, spark region, and spark The refracted angle 1S associated with the refractive index 
power A -spark” provides a fixed amount of energy 
(spark energy) in a region (spark region) centered around 
the focal point for a fixed amount of time (spark time). 

When the strength of a spark reaches a threshold, or the 
threshold breakdown intensity is satisfied, a laser 
supported combustion (LSC) or detonation (LSD) is 

initiated. , 

That means enough seed electrons are produced 

through the likes of associative ionization reactions sue 
that the subsequent optical breakdown (ignition) is 
possible. It is found that an ignition is sustainable when 


A ray may change its propagating direction due to the 
inhomogeneous refractive index within the hot plasma. 
The index of refraction is taken from Edwards and 
Pleck: 18 



1 - 


n = 




ni'C-n ]_ 


through the Snell’s law: 

rtj sin 0; = sin $2 

where m and n; represent the refractive indexes tor two 
different control volumes, and 0i and 0; are the inu ent 
and refracted directions with respect to the norma 
direction of the interface between two volumes. 
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experiments. The CO; laser delivered up to 400 J single- 
pulses at a pulse width of 18 ps. The laser pulse energy 
was measured with a calorimeter and the uncertainty o 
the delivered energy was estimated to be ± 10 J. Severn 
variations of the basic Laser Lightcraft design (Model 
200 series) similar to those described in Ret 8 were 
examined and the test results of the 6061-T6 < all 
aluminum Model #200-3/4 vehicle and Model A vehld ® 
are chosen for this study. The vehicle shape of Model A 
can be found in Ref. 7. The impulse measurements were 
conducted with a pendulum apparatus. This technique 
employed a velocimeter coil, which was used in previous 
work to determine the impulse imparted to a flat plate 
using the PHAROS III laser at Naval Research 
Laboratory. 30 The Laser Lightcraft was suspended and 
suitably weighed before being subjected to a single pulse 
of energy from the PLVTS laser. The uncertainty in the 
impulse measurement is estimated to be 1% or better 
Details of the laser, measurement technique, and 
pendulum apparatus can be found in Rets. 7 and 8. 

f nmpiitat ipnnl Grid CiPnflratlQII 

Figure 1 shows the layout of a computational grid (for 
Model #200-3/4). As described in Refs. 7 and 8. the 
Laser Lightcraft vehicle consists of a forebody (nose), an 
annular shroud, and a PDW engine (parabolic optic L 
Only half of the grid shown in Fig. 1 is actually solve 
due to the axisymmetric formulation. An eight-zone 
structured grid was generated first using a grid generation 
code called UMESH. 3 ' The grid is then fed to an 
auxil.arv program such that an one-zone unstructured- 
arid is constructed for the actual computation 
Quadrilateral elements are used such that the wall 
boundary layer can be computed it needed. Hig gri 
density is used in the inner shroud region for capturing 
the optical breakdown and pulsed detonation wave 
processes. 


Poc.ilK and Discussion 

A series of computations have been performed tor 
average laser pulse energies of 75, 100. 150, 200, 
and 400 J for Model #200-3/4, and 400, 600. and 800 
for Model A. Maximum computed temperatures ranging 
from 51 000 to 60.000 K fall within the computed post- 
shock temperatures of 20.000 to 60.000 K using Park s 
multi-temperature air chemistry. LSD occurs in all 
cases with maximum Mach number reaching -.3. LSL 
then ensues as the plasma front expands, and the Mach 
number decreases to subsonic value. Figure _ shows the 
computed temperature contours and laser beam traces at 
an elapsed time of 10 ps and an average laser pulse 
energy of 400 J. It can be seen that the laser beam 
reflects specularly on the optical surface and focuses onto 
a focal “point” on the shroud where the breakdown ot air 
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starts. Ten laser rays are plotted for clarity, whereas two 
hundred rays are actually used in the computation. 

It is probablv best to recount the phenomenon ot an 
optical detonation wave, which is a detonation wave 
oenerated by the optical breakdown of a medium with an 
incident laser beam, at this moment. This was put 
forward by Raizer 1 ' in detail in his book. Essentially, a 
strong shock wave is generated in the region where t e 
laser radiation that produces the plasma is absorbed and 
heat is deposited very rapidly. This shock wave travels 
through the medium, heating and ionizing it such that the 
medium becomes capable of absorbing more laser 
radiation. The laser energy is deposited in successive 
layers of a medium adjoining the front of the shock wave 
subjected to the laser radiation. These layers themselves 
thus become energy sources that maintain the shock 
wave. Hence, the shock wave moves along the optical 
channel in the opposite direction to the laser beam. Thus, 
during the initial stages the shock wave is maintained by 

the laser beam and does not decay. , 

In 2 the optical breakdown is being ted by the 
laser ene°rgy and the plasma front grows. Notice the rays 
are allowed to bend as the index of retraction varies with 
the expanding plasma. As the local electron density ot a 
plasma front reaches the critical value such that plasma 
resonance occurs, the local plasma becomes opaque and 
the regular laser energy absorption is greater on other 
points of the surface of the plasma front. This process is 
nonlinear and may explain the irregular shape of t e 
plasma front. It is anticipated that eventually the laser 
beam is totally reflected as the plasma resonance 
phenomenon spreads over the entire plasma front. The 
phenomenon of plasma resonance disappears as the 
plasma front expands and the electron density decreases. 
The computed optical breakdown phenomenon appears 
to a«ree with the description given by Raizer. Notice 
that °the spreading of the shock wave outside the optical 
channel results in gradual attenuation because there is no 
energy to maintain it. Thus, as the shock wave expands 
in the direction opposite to the laser beam, it is weakly 
maintained by laser radiation when it is not in a state ot 
plasma resonance. 

Figure 3 shows the temperature contours at an elapsed 
time of 20 us. At this time, the laser beam is turned oft 
since the pulse width is fixed at IS ns. thus the shock 
wave is not being maintained by the laser radiation. The 
computational results show that a pulsed detonation wave 
usually takes place at approximately 10 to 20 jis after e 
start command. The higher the pulse energy, the earlier 
the generation of a pulsed detonation wave. Part ot this 
propagating detonation wave will hit the parabolic optic 
and reflect off the surface. The reflected and non- 
re fleeted portions would combine to make a propagating 
Mach stem. It should be noted that the phenomenon of a 
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propagating Mach stem on the optical surface is very 
s Jlar to that of a propagating 1 Mach^sc ^beU 
nozzle during the start-up transient. This is becaus 
the aerospike optic geometry is an "inside-out version 

that of a contoured bell nozzle. 

As the Mach stem moves past the end of the shroud, 
the part that is attached to the inner-shroud starts to wrap 
around the outer-shroud. As the wrap around progresses 
the shock wave weakens and grows like an ellipsoi . 
forebody attaching “leg” eventually develops as shown 
in the pressure contours of Fig. 4 at 1050 ps. Eventually, 
the "leg” leaves the nose and the ellipsoid expands to 
infinity. That completes the propulsion physics of the 
Laser Lightcraft for a single pulse of laser energy. 

The coupling coefficient is the generated impulse (N-s) 
divided by the delivered average laser single-pulse 
energy (MJ). It is a performance measure unique to the 
pulsed laser propulsion, and conceptually similar to the 
specific impulse (generated thrust divided by propellant 
weight flow rate) of a chemical rocket. Figure 5 shows a 
comparison of model predicted and experimental 
measured coupling coefficients for Model #200-3/4 The 
model with non-equilibrium radiation produced slightly 
better comparison than that without non-equilibrium 
radiation. The scattering of the measured data jis 
attributed to the measurement of the laser energy, w ic 
varied from shot to shot (± 10 J) and an average over 
many pulses was used. The model (without non 
equilibrium radiation) predicted coupling coefficients 
a«ree reasonably well with those measured in the 
laboratory. Figure 6 shows a comparison of model 
predicted and measured impulses tor Model A. The 
model predicted slightly higher impulses at laser pulse 
energies from 600 to 800 J. It is anticipated that the 
comparison improves with inclusion of non-equilibrium 

calculation. , . „ 

Fmure 7 shows the predicted maximum electron 

density (N«.) as a function of elapsed time for Model 
#200-3/4. It rises sharply within about 5 ps, reaching its 
peak, then drops with time. The earliest rise and t e 
highest peak value occur for the highest laser pulse 
encruy (400 J). as expected. The computation indicates 
the occurrence of plasma resonance for all cases since t e 
maximum electron densities exceed the critical value tor 

about 25 to 45 ps. 
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detonation wave propagation agree well with those 
... I rnrxHM nrefilCted COUDline 


detonation wave t . . .. 

described in the literature. The model predicted coup i g 
coefficient for a Model #200-3/4 Laser Lightcraft and the 
predicted impulse for a Model A Laser Lightcraft agreed 
reasonably well with those measured. 
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fnnrlusions 

A computational plasma aerodynamics model has been 
developed to study the propulsion physics of an 
experimental Laser Lightcraft. The model development 
.s based on a building block approach such that the 
model can be improved continuously with improved 
understanding of the physics. The model predicted user 
propulsion physics such as the optical breakdown an 
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